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Abstract 


An elementary derivation and a complete description are given of an 
algorithm for Interpolation over a plane triangle when function values 
and first partial derivatives are given at the vertices. The method 
gives continuity witl. neighboring triangles. 

The interpolation method is mathematically equivalent to one that has 
been discussed previously in the literature; however, the algorithmic 
form given here is more efficient than has previously been described. 
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Compatible Interpolation Over A Triangle 


1 . Introduction 


The problem treated in this report has been treated by numerous 
authors. See Birkhoff and Mansfield [1974] for extensive discussion of 
this and closely related problems and for other references. inter- 

polation over triangular grids has application in structural analysis via 
finite element methods and in the computerized representation of surfaces 
for computer aided design. 

The problem also arises in diverse scientific and engineering fields 
where it is useful to be able to construct a smooth surface that passes 
through a finite set of observed or computed values of some function 
z = f(x, y). In these data-fitting applications, the desired end-product is 
often a contour plot of the interpolated function. 

In this report, we give an elementary derivation and a complete algor- 
ithmic description of an interpolation method that is mathematically equiv- 
alent to one that is mentioned by Birkhoff and Mansfield, [1974], and spec- 
ified in detail by J. -J. Goel, [1968], Goel attributes the method to Clough 
and Tocher [1965] and Zienkiewicz [1967], 

The algorithmic form in which the method is given here is more 
efficient for solving the interpolation problem than the form given by Goel 
[1968]. It is my present conjecture that one cannot expect to discover an 
algorithm for this problem that is significantly more efficient than the one 
given here. Other reports, yet to be written, will deal with the integra- 
tion of this interpolation algorithm into a set of subprograms for constructing 
a triangular grid, Lawson [1972], and then doing look-up, C^ interpolation, 
and contour plotting for a function, z = f(x, y) , whose values are given at a 
finite set of points. 

The author thanks Dr. Fred T. Krogh for numerous fruitful dis- 
cussions during the exploratory phase that preceded the writing of this 
paper and for a critical reading of the paper that produced numerous im- 
provements. 
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2 . 


The problem of C ~ compatible interpolation over a triangle 

Assume that values of a function, f, and its first partial derivatives, 
f^ and fy, are given at the three vertices of a triangle, T, in the (xy)-plane. 
We wish to define a function w(x, y) for (x, y) in the triangle, T, that will 



With nine items of data being given we may anticipate that the inter- 
polation method can be required to be exact for all polynomial functions 
of degree up to two, but not for all cubic functions, since the set of quad- 
ratic functions in two variables is a six-parameter family while the set of 
all cubics is a ten-parameter family. We will impose the requirement 
that the method be exact for quadratic polynomial data. 

Furthermore we want the interpolation method to have the property 
that if it is applied to two adjacent triangles having an edge in common 
then function values and the first partial derivatives of the two interpolated 
functions will be identical along the common edge. Thus the method can 
be used for interpolation over a triangular grid, and the surface defined by the 
totality of the locally interpolated functions will have continuity over 
the entire region covered by the grid. 

A convenient way of assuring that the interpolated functions on adjacent 
triangles have the same values and first partial derivatives along the com- 
mon edge is to require that the values and first partial derivatives of the 
interpolated function along any edge must be determined only by the data 
given on that edge, i. e. the data given at the vertices at the ends of that 
edge. 
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Relative to a particular edge, say Sj, the given partial derivatives at 
vertices V 2 a^nd can be rotated to give partial derivatives tangential to 
the edge and normal to the edge at W ^ Vj, A fairly natural approach is 
to define values of w along side Sj by Hermite cubic interpolation matching 
the required function values and tangential first partial derivative values 
at V 2 «^nd Vj.and to define the first partial derivative of w normal to side 
Sj to be a linear function along side S^, taking the required values at 
and Vy 
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Some observations on one ^ dimensional Hermite cubic interpolation 



Consider Hermite cubic interpolation on the unite interval 0 ^ x ^ 1, 
with given data £q and at x = 0 and £j and £^ at x = 1. The cardinal 
£unctions £or this data are 






The interpolated cubic polynomial is given by 


( 5 ) w(x) = fo<Po(x) + £JPq(x) + f cpi(x) + ^^(x) 

Note that all four of the cardinal functions in this formulation are cubic 
polynomials. In our interpolation problem over a triangle we shall find 
that cardinal functions of degree higher than two introduce relatively 
large increases in computational complexity. Thus it is useful to note 
that the solution to the one -dimensional Hermite cubic interpolation 
problem can be rearranged (in various ways) to involve at most one 
cubic basis function. 

For example, we can use the four functions: 


( 6 ) tj(x) = 1-x 


( 7 ) tjCx) = X 


6 
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( 8 ) ^ 2 ^*) = 


and 


( 9 ) tjCx) = 2x(x~ j-)(x-l) 
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The same cubic polynomial as is defined by £q (5) can then be constructed 
using the formulas 


( 10 ) 


m = fi - fo 


£1 - f« 


(11) w(x) = fQi)j(x) + fjifj(x/ + 2 — 1 ) 2 ^ 


f' + f' - 2m 


K(x) 


Eq (11) is easily derived by constructing the formula in two stages. 

The first two terms clearly provide the linear interpolant that matches the 
data Iq and f^ at 0 and 1 respectively. This linear function has a 
slope of m = f j - f^. 

Thus, after subtracting this linear function, the remaining problem 
is to determine a cubic function having zero values at the endpoints and 
slopes of fQ - m at 0 and f'^ - m at 1. Since 1)2 slopes of 1 and 
-1 at 0 and 1 respectively and has a slope of 1 at both 0 and 1 it 
follows that the function 

t2(x) +-2 ^ *3(x) 

will fit the residual data. 

Combining this two-stage procedure into a single formula gives Eq (11). 

An analagouh' approach of combining a linear interpolant with a quadratic 
part and then a cubic part will be de'-cribed in Section 5 for the i'lterpolation 
problem over a triangle. Before commencing this derivation, however, we 
must introduce definitions and notation for the coordinate systems we will 
use over a triangle. This is the subject of Section 4. 
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Wr^': 


4. Coordinate systems over plane triangles 


Let T be a plane triangle having vertices / ^ 2 ' ^3* 

convenience, we will thir ’< of these vertices as being labeled in counterclock- 
wise order. Final results are not sensitive to this assumed ordering. 

The indices in the formulas to follow take the values 1, 2, and 3. 
Index arithmetic is to be interpreted as being cyclic over these three values. 
For example, if i = 3, then i + 1 = 1 and i + 2 - 2. 

Let S. denote the side of the triangle opposite to vertex V^. Let 
e denote the interior angle at vertex V., measured from side S.^2 

4 . 



In a Euclidean (x, y)-plane let (Xj, y^,) be the coordinates of vertex 
V., i = 1,2,3. 


Introduce directed edge vectors e^^ with components u^^ and v^ defined 


by 


(13) = 


m m 


m 

u. 

1 


’‘i+2 ■ ’‘i+1 

V. 

1 


n+2 - Vi+l 

m m 

1 

- 


; i 
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Denote &e Euclidean length of side by 
t - ii_ II _ #.2 . 2.1/2 


A = ¥i\\ = («i + V. ) 


i = 1,2,3 


Let c. denote the inner product (dot product) of the two edge vectors directed 
away from vertex V^, i. e. the inner product of and • 


(14) c, = Dot(VjV,^.,. VjVj.j) = Dot(e..,, 

= -‘"i-i”!*! * vri+i* 

— COS I “ 1, 2, 3 

We note in passing that the c.'s and /^'s are related by the equations 




1 = 1.2, J 


-l^ - c - c 

^i+1 Vl i-l % 


i = 1,2,3 


2 c. = i... + if , - /f 
1 x +1 1-1 1 


i = 1,2,3 


Let 6 denote twice the (signed) area of T. This quantity is repre- 
sentable as the scalar cross product of any two edge vectors directed away 
from a common vertex. 


10 
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(18) 6 = Cross (VjV.^j, ^ ^Vl’ "®i+l) 


= Cross (e.^j, e._j) 


= det 


Vl Vl 


V.,, V. , 
1+1 1-1 


= /. sift 0. 

1-1 i+l 1 


ft-j.iV. , - V. , ,u. . 

1+1 1-1 1+1 1-1 


i = 1,2,3 


Tangential and Normal Coordinates 


Relative to side of the triangle, we introduce an orthogonal 
coordinate system using coordinates (t., n^) where t. is measured along 
side from and n^ is measured positive in the inward normal direction* 



The variables (t^, n^) are related to the variables (x, y) by the 


equations 


u. V. 

1 1 


-V. u. 

1 1 


y-n+i 


i = 1.2, 3 


and by the * tverse equations 
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1 




u. -V. 

1 1 


V. u. 

1 1 


i = 1.2.3 


From £q. (20) we obtain the partial derivatives 


ax/atj = Uj/ij 


ax/an. = -V.//. 

1 11 

dy/at. = v./y. 

ay/an^ = “j/Zi 


i = 1, 2. 3 


Barycentric or Areal Coordinates 

Let P be a point with coordinates (x, y). Define the three barycentric 
(or areal) coordinates of P by: 


= Cross V.„P) / J 

= Cross (Cj, Vj_j.jP) / 6 



U. 

} 

*‘*j+l 

det 




v. 

J 


[«j(y 

■’'l+l' ■ 

v.(x-x.^l 


= 6 [u.(y-y.^j) - v.(x-x.^j)J j = 1,2.3 

Note that the quantity computed by the cross product in the formula 

for r. is twice the area of the triangle formed by side j and the point P. Thus 
J 

the sum of the three cross products used to compute rj, sind r^ must be 
twice the area of T. Therefore, with the normalization factor ( ^ appearing 
in the formulas it follows that 


12 
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(23) 


1 


-i 


( . 


■1-- 


Ti + T2 + rj = 


The barycentric coordinates are the unique set of numbers having 
unit sum and representing P as a linear combination of VjV 2 * sind thus 


r 1 

X 

3 

X. 

J 


M II 
II 

y. 

/j 


Each barycentric coordinate r^ is the unique linear function of (x, y) 
that is zero along the line determined by side and takes the value 1 
at the vertex V., 

For points inside the triangle T we have all r ^ a 0 while for points 

outside there will be some r. < 0 . The barycentric coordinates of the 

J 

vertices are: 

Vj ~ (1,0,0) 

V 2 (0. 1,0) 

V 3 ~ ( 0 , 0 . 1 ) 

Using Eq. ( 22 ) we may compute partial derivatives as follows: 

(24) 


ar./ax = -V./6 

J J 


br./by = u./i 


j = 1,2.3 


Using Eq. (21) and (24) we obtain further partial derivatives; 
(25) sr^/at. = (arj/ax)(ax/at.) + (brj/ay)Oy/atj) 

= + UjV.) / (/j6) 

= Cross (Cj, e^) / (/^ 6 ) 
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= Jl//. i£i-j = l 

(- 1 // 

(26) = (arj/ax)(^x/an.) + (arj/ay)(ay/anj) 

= (VjVj + UjUj) / (/j6) 

= Dot (e., e.) / 

tl 

/f/(ii6) ifi = j 

-Ci^.i/(^i6) if i-j = 1 

For more convenient reference, we organize the results of Eq. (25) and (26) 
into tables as follows: 

Table 1 . Values of 3r./at. 

J 

arj/ 

arj/ 


/atj /at2 /atj 


0 

1//2 



0 



-1//2 

0 



Table 2 . 


Values of afi/an- 
















r 




Constructing a solution 


First convert the given partial derivative data at each vertex to partial 
derivatives in the directions of the edges meeting at the vertex. For i = 1. 2, 
and 3, let h^ and respectively denote the values of the first partial deriva* 
tive with respect to t^ (the tangential direction along side S^) at t^ = 0 and 


tj = I (i. e. , at the vertices and 


'2X h 




The values of h^ and are given by 


*'i = <“i ‘x, i+1 + n 'y, 
''i = <"i ‘x, 1-1+^ «y. 


i = 1. 2. 3 


Henceforth, we may regard h^, and k^ for i = 1, 2, 3, as defining the data 
to be interpolated. We proceed by analogy with the discussion of the one- 
dimensional problem in Section 3 . A linear function interpolating the func- 

tion values f^, i = 1, 2, 3, is given by 

o,<‘> = ri<i + 


Along side S.. this function has slope 

(2») m. = (fj., - £.^j) /«. 
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with respect to the tangential variable t. . Subtracting the vertex values and 

( 1 ) ^ 

partial derivatives of («' ' from the given data, we are left with the residual 
problem of interpolating vertex values of zero and tangential partial deriva- 
tive values of h.-m. and k.-m., i = 1, 2, 3. 

IX 11 


Let us temporarily restrict attention to one side, say side Sj. On side 
Sj we have r^ =0, and tj = Note that the quadratic function 

reduces to r 2 (l~r 2 ), and the cubic function r 2 ^^ r reduces to 
2r2(r2“"2 )(r 2 "l) along side Sj. [Compare with 1^2 +3 Section 3 . ] 

It is easily verified that the partial derivative with respect to tj of the scaled 
function Jtjr 2 rj has the values 1 and -1 at the vertices V 2 and respectively. 
Similarly d[t jr 2 r 2 (r 2 "r 2 ) ]/atj has the value 1 at V 2 and V^. 


Thus the function 


(2) ^r*^l . . 

(29) it>\ = — 2 — ^1*^2^3 5 ^ir,r,(r,-r,) 


'1*2‘3'‘2 ‘3' 


( 2 ) 

satisfies the residual interpolation requirements on side S^*, i.e. , u>'^ ' 

has zero values at V 2 and and its partial derivative with respect to tj has 
the values h^-m^, and kj-m^ at V 2 and respectively. 

( 2 ) 

On side S, we have r, = 0, and thus (ju\ ' and its tangential partial deriva- 

^ ^ ^ (2) 
tive are zero there. Similarly, since r^ = 0 on side S^* oi)j and its tangential 

partial derivative are also zero on side S^. 

( 2 ) ( 2 ) 

By appropriate cycling of indice s^ define f mctions u)i ' and u)' ' analagous 

(2) . ^3 

U)j 


to 


h.-k. h +k.-2m. 

( 30) Jt.r. . .r. , + ^A.r. ,r. ,(r. , - r. ,) 

' ' 1 2 1 1+1 1-1 2 1 1+1 1-1' 1+1 1-1' 


for i =1, 2, 3 


16 
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It follows that the function 


( 31 ) 


(U 



3 

E 

i=l 



interpolates the nine items of data f., h., k., i=l, Z, 3. 

Note further that u) is exact for quadratic functions since if the data 
f^, h^i i=lt 2, 3, arise from a quadratic function it will follow that 

h.-m. = -(k.-m.) i - 1, 2, 3 

1 1 ' 1 i' 

so the coefficients of the cubic terms in ud vanish leaving o) the unique 
quadratic function matching the given data. 

The function o) has a defect however. We require that the partial 
derivative normal to any side must be a linear function along that side. The 
cubic functions ^ i-l^*^i+l'*^i-l^’ have this property and 

thus, in general, ^ does not. 


The remedy, described in Goel [196^, is to introduce correction 
functions p., i=l, 2, 3. The function is required to be zero on all three 
sides of the triangle. It follows that its firsv partial derivatives in all direc- 
tions at each vertex are zero. It is further required that the normal derivatives 
of relative to sides and be zero on those sides respectively, 

while the normal derivative of p^^ relative to side is to be a quadratic 
function along that side. Specifically we can require that 


(32) 




on side S 


i 


By adding appropriate multiples of p^, p^ to a cubic function, 

such as i^ 2'3^'2''3^' construct a function whose normal partial der- 

ivatives on each side are linear functions along the respective sides. 
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This point will be further developed in the next three sections. In 
Section 6 we determine the multiples of p^, p^, and needed to correct the 
function Sections 7 and 8, we discuss two distinct sets of 

functions having the properties required of the p.'s. 
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6 . 


Correcting the normal derivatives of gj = r 2 r 2 (r 2 "rj) 


Define 


8l = 

The partial derivatives of gj with respect to the r/s are 

9gj/Srj = 0 

Sgj/dr2 = 2r2fj - r^ 

Using the expressions for dr./9n. given in Table 2 , and evaluating on 

^ J 

the indicated sides we obtain: 


(33) 


[3(C3-C2)r3 + 2(2c2-C3)r3-C2]/(iiM 


on side S, 


(34) &gj/9n2 = " /t 

on side $2 

(35) agi/Sn3 = ljr^/6 

on side S, 


We assume the availablity of functions p., i=l» 2, 3, which 
have zero values along all edges and satisfy 

( “ J '* 

SPi/an, I = ( 


ap./3n. = < 

^ ^ -J C lo ifj=i-l or j=i+l 

on side S. 

1 


We wish to determine coefficient ck^ji 0fj2» ®13 swch that the function 


t \ - «1 * 11^1 * 12^2 ^ “ l 3^3 
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will have the property that for j = 1, 2, 3, the normal derivative 9gj/9nj 
is a linear function along side Sj. Comparing the quadratic terms in 
Eqs. (33), (34), and (35) with the quadratic term in Eq. (36) it follows that 
the appropriate values for the ®jj's are 


a 


11 


3(C3-C2)/Jti = 


and 


“12 = ■* 


“l3 " ' 


Using these values of the obtain the following equations 

showing the linear character of the normal derivatives of on the respective 

sides: 


a(Aii'i)/anj 


” t(c2+C2)*’2"C23/4 


Ion side S 


1 


= (C3r3-C2r2)/* 

'3 . '2 

tanST tand. 




- “ A j f 2( “ ”^1^2^3^^ 


on side S. 


•r^/sinSj 




on side 


= = r2/sin02 


Collecting the results of this section and cycling the indices ap- 
propriately we define 


20 


JPL Technical Memorandum 33-770 







(37) 


= ^i+iVi<vrvi> + 


^f-i> 


Pi"pi+i^Pi-i 


i = 1. 2, 3 


Each function g. has the same interpolatory properties as the simpler cubic 

function j) three vertices and has the additional 

property that /dn. is a linear function along side S. for all i and j. 

^ J 3 

Therefore for our complete interpolation formula we replace 
Eq ( 31) by 


3 r h.-k. 

r = E Ir.f. + — 5 — ^ i.r. ,r. 

I 1 1 2 X x-1 iH 


h.+k.-2m. 
11 1 


- ‘a1 
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7. 


A set of rational correction functions, 


Define 


(39) 




1 . 


This set of rational functions is discussed by Goel [1968] 
use in this context to Zienkiewic% [19S7]. 


2. 3 

who attributes their 


For convenience consider the single function 
Pi = I 

Over the triangle, T, the denominator of p, vanishes only at the vertex V 2 » 
where r 2 = 1. and at V^, where r^ = 1. These are removable singularities 
however since, for example, at V 2 the numerator has a third order zero 
(rj - 0 and r^ = 0) and thus Pj has a second order zero at V 2 . Similarly 
Pj vanishes to second order at V^. 

At all other edge points it is clear that p^ vanishes because it contains 
Tj, r 2 » and r^ as factors. 

To verify that the normal partial derivatives of pj have the necessary 
properties compute 


= (2-r2)Pj/Cr2(l-r2)3 

apj/Srj = (2-r3)Pj/[r3(l-r3)] 

Then using expressions for dr^/dnj from Table 2 


, one finds 
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i 


i 


> 



dp j/arij 


= r3(l-r3)ij/ft 


on side S 


1 


dp j/dn2 


= 0 


on side S> 


9pl/9n3| 

ion side 


It is thus verified that the rational functions of Eq(39) can be used 

as the coriection functions p. described in Section 5. 

1 


To compute values of p^^, i = 1, 2, Z, assume values oi rj^, 
i = 1, 2, 3, and 


(40 ) 


®. = r..,r. , 
1+1 1-1 


i = 1. 2, 3 


are given. One then computes 


♦ := rj9j 

:= 1-r. i = 1, 2, 3 

(if any = 0 branch to handle the special trivial case of 
interpolation at a vertex} 


P 


i 




i+l 



i = 1,2,3 


Thus the computation of the rational p^'s requires 7 multiplications, 3 additions, 
3 divisions, and 3 zero tests. 
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A set cf piecewise cubic correction functions, 


Define 


' *‘i^^*‘i+l*^i-l^ *‘i ” ““ *^2’ ^3^ 


(<>) 0 , 


■’i+l'-'l+l + 




if r^^j = mintrj, r^, r^} 
if r^_i = min(rj, r2» r^} 


for i = 1, 2, 3 


Each function p. consists of a set of three cubic functions which match with 
1 1 

C continuity along internal boundary lines connecting the vertices to the 
centroid of the triangle T. This set of functions is discussed by Goel [1968] 
who attributes their use in this context to Clough and Tocher [1965], 

Along side Sj we have r^ = min{rj, r 2 » r^} and thus 

?1 = rj[6r2rj + r^lSr^ - 3)1/6 

In this region the derivatives 9pj/^r^ are given by 

apj/drj = [6r2r2 + IStj - 6rj]/6 

aPl/ar2 = rjr3 


9Pl/5r3 = r^r2 


Using the expressions for dr^/dn^ from Table 2 > we find 
Spj/dni j - 

ion side 


It can also be verified tl,''*' 


24 
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9p j /Sn2 


and 


= 0 

on side 


apj/3nj j 

Ion side 


0 


Thus the piecewise cubic functions of Eq (41) have the properties need for 

use as the correction functions p. of Section 5. 

1 


If Eq (41) is used in a computer program the usual situation will be 
that for one set of values of r,, r 2 » and r^ the program must compute values 
of p^, P 2 > a^nd p^. In this context the following restatement of Eq(4l) is 
useful: 


Let m be an index such that r = min(r,, r^, r,}. Then 
... m i c i 

p = r [6r ,,r ,+ r (5r -3)]/6 

^m m*- m+1 m-1 m m ' 

P^., = r^(-r +3r ,)/6 

m+l m m m-1 

= + 3 r ,)/6 

p , m m m+l 

^m- 1 

To compute values of p^, i=l, 2, 3, given r^, i=l, 2, 3, and 


i ■ ’^i+l'i-l 


i = 1, 2. 3 


the following steps can be used: 


Find m such that r =minfr,, r.,, r_} 

m ICO 


a 

b 


1 

2 

\ 

3 


r 


2 

m 


r 


m 
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m 

:= r (® + 4 a) - a 

m^m 3 

^m+l 

■= 

*^m-l 

• *>> 


Two compares are required to determine m. Counting these as additions, 
the computation of the piecewise cubic correction functions, p., requires 
seven multiplications and six additions. 


26 
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9. Transforming formulas to algorithms 

Assume that data (x., y., f., f f •)» i - 2, 3 is given, 

X X 1 XyX ^yX 

where (x., y^) specifies the coordinates of the vertex and (f., f^ f^ 
specifies the value of the function and its first partial derivatives with respect 
to x and y at the vertex V.. The vertices should not be colinear. Either 
counterclockwise or clockwise ordering of the vertices is acceptable. 

Further assume a coordinate pair, (x, y), is given at which an inter- 
polated value, w, is to be computed. The point (x, y) should not be exterior 
to the triangle T having vertices Vj, V 

We will describe the interpolation algorithm in three phases. Phase 1 

2-1 

will compute the preliminary quantities, u^^, v^, jl^, 6 , r^^, and cp., i = 1, 2, 3. 

Phase 2 will be the computation of the p^'s using either the method of Section 7 
or of Section 8 . Phase 3 will complete the Interpolation. We will describe 
three versions of Phase 3. 

All index expressions such as i+1 and i-1 must be in,erpreted 
cyclicly so that the resulting index value is 1, 2, or 3. The name "det” is 
used to indicate computation of the determinant of a matrix. Phase 1 of the 
algorithm proceeds as follows: 

Begin Phase 1 


u. := 

X 

Vl ■ 

*i+l 

V. : = 

X 

Vl ■ 

n+1 


2^2 
U. + V. 
X X 

det 

'^1 "2 



'^l '^2 



{Test for the error condition, 6=0, indicating colinearity of the 
of the vertices. } 
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X = x-x. 


y '= 


*^2 


y-Yi 


6 ^ det 


u- 


a 


• 


u. 


= i det 


X 

/W 

y 


rj := l-(r2+ r^) 


{If one wishes to test for the possibility that (x, y) is exterior to the 
triangle the test can be made here. The condition ^ ^or any i 

indicates that (x, y) is exterior to the triangle. ] 

End Phase 1 

Note that Phase 1 requires 17 multiplications, 16 additions, and one division. 

Phase 2 consists of the computation of the p/s using either the 
method v" Section 7 or of Section 8. We proceed to the discussion of Phase 3. 

The divisor in Eq. (27), and (28) will cancel with the multiplier i. 
in Eq (38). Thus instead of computing h. and k. we will compute quantities 
h. and k. such that h^ = ^i^^i ^i ~ 

Version 1 of Phase 3 is the following: 
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Phase 3, Version 1 


'll. 


k. 

1 


/%/ 

«i 


:= u. f . , , + V. f . , , 
1 X, 1+1 1 y, 1+1 


;= u. f . , + V. f . , 
1 X, 1-1 1 y, 1-1 


2 2 
Vi - ^i-i 


:= + ^ 12 


Pi ■ Pi+1 + Pi-lj 


i = 1, 2. 3 


w : = 


1 yV' ^ 

* #■« 1 


1 


Z f.r. + - 5 - (h. - k.) {p. + C-|<h. + k.) - £. , + £.^, ] %. 

i 1 2 ' 1 i' *'2' 1 x' 1-1 1+1-* ®i 


Phase 3, Version 1, requires 36 mulUplications, 42 additions, and 3 divisions. 
This includes six multiplications by one hal£. The computation can be re- 
arranged so that there is only one multiplication by one hal£. This leads to 
what we will call Phase 3, Version 2 which requires 31 multiplications, 43 
additions, and 3 divisions. 


which 


Version 2 di££ers £rom Version 1 only in the expression £or w 
changes to: 



From this point it takes only a little more rearran'*'^ nent to obtain a 
£ormulation which explicitly uses cardinal £unctions £or the give i data 
(£., £ ., £ .), i= 1, 2, 3 This £ormulation, which we will call Phase 3, 

t X, 1 Yi 1 

Version 3, uses the same number o£ multiplications, additions, and divisions 
as Version 2. 
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R- •' 


Phase 3, V e rsion 3 

^ '^^2 ™ Pi'f*i+1 ^ Pi-1 

Ji» 

1 


Pi := gitfi 




g- ■ 9- 


i = 1. 2, 3 


“i '= 'i + gi-1 • 8i+l 

h ’= ViPi-1 ViVi 


'^i-iPi-i + \+iVi 


ill 'i“i " 2 1 "i " V. ‘ "i ) 


i = 1. 2. 3 


Version 3 is particularly efficient for the case in which there are a 
number of functions to be interpolated at the same interpolation point, (x, y). 
In such a case only the final formula for w must be recomputed for each 
function to be interpolated. 


The explicit computation of the cardinal functions, a-i -^3^, and 
•| V- , as provided by Version 3, is needed in applications in which the quantities 
f , f and f . are unknowns to be solved for. This is the situation in the 

*1 Xf X yi 1 

finite element methods for solving partial differential equations and in the fitting 
of a smooth surface to noisy data. 
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10 . 


The method described in Goal C 19681 


For comparison with the results of Section 9 we will give a brief 
description of the interpolation method given in GoSl [1968]. This method 
is attributed by Go81 to Clough and Tocher [1965] and Zienkiewicz [1967]. 

We refer primarily to Eq (28) and (34) of GoSl [1968], The fol- 
lowing change of notation will convert GoSls symbols to ours. 


GoSl's notation 

Notation of this paper 

x' 

*•2 

y' 

'3 

1 -x'-y' 

^1 

"i 

-Pi/6 

“i' fi’ ''i 

Ofj, Y; 

^ 

“i* Pi* ^i 

*i’ 

A 

6/2 

Li. C. 

^i’ ‘'i 


We will indicate precomputation of the common ^iibexpressions in 
GoSl's formulas in order to provide a basis for obtaining a realistic operation 
count. The quantities computed in Phases 1 and 2 [see Section 9 ] are all 
needed for GoSl's formulas so we will assume these computations have been 
done and proceed to describe Phase 3. 
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I 


Phase 3, GoSl 


\ := ^(2r.,,-3) 

^ := Pi/^i 

^ == 'i«Pi 

1 + + Vl - 2 [Pi ■ 2 (Pi+i + Pi_i - f)3 

'iVl ■ Pi" ^Pi+1 s i = 1. 2. 3 

iVl ~ Pi"^ ^Pi+1 ■ ^Pi-1^ 

«i == “i " <Pi+i + Pi-i3<ViVi + ViVi> 


/\ 

“i = 


A 

P; : 


A 

Yi ^ 


= r. 


Pi 

^ A. 1 , . ^ 

< 

1 

1 

Yi 

A A 1 , - 

+ u. ,p. , 
1- 1 1- 1 

w 

»i * *x. i * V i "i’ 

1=1 



Assuming the quantities •I’e, 3p^, and 5p^, would be computed only 
once eacht the operation count for Phase 3, Gofil, is 71 multiplications, 

81 additions, and 3 divisions. 


It can be verified, by the appropriate tedious algebra, that a^, 
and V- S’® defined above are identical with snd ss defined in 

Phase 3, Version 3 [Section 9 ]. Thus the interpolated value, w. is the same 
by either method. 
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Summarv of ODeration counts 


Phase i 


Phase 2 using rational pJs 
Phase 2 using piecewise 


cubic p^'s 


Phase 3, Version 1 
Phase 3, Versions 2 or 3 
Phase 3, GoSl 

Totals using piecewise cubic 
Pi's : 

Version 1 
Versions 2 or 3 
Goel 


Multiplications 

Additions 

Divisions 

17 

16 

1 

7 

3 

3 

7 

6 

0 

36 

42 

3 

31 

43 

3 

71 

81 

3 

60 

64 

4 

55 

65 

4 

95 

103 

4 


If we weight the multiplications, additions, and divisions in the ratio 
2:1:6, then the operation count for Versions 2 or 3 is 63% of the count for GoSl's 
version. 
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